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We propose a new method for calculating total energies of systems of interacting electrons, which 
requires little more computational resources than standard density-functional theories. The total 
energy is calculated within the framework of many-body perturbation theory by using an efficient 
model of the self-energy, that nevertheless retains the main features of the exact operator. The 
method shows promising performance when tested against quantum Monte Carlo results for the 
linear response of the homogeneous electron gas and structural properties of bulk silicon. 
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PACS numbers: 71.15.Nc, 71.45. Gm, 71. 15. Mb 

Density-functional theory (DFT) Q is a powerful 
method to calculate ground-state properties of electronic 
systems. In its standard Kohn-Sham form Q, the sys- 
tem of interacting electrons is mapped onto a system 
of non-interacting electrons moving in an effective local 
potential. This potential, however, exhibits some non- 
analyticities in its dependence on the electron density, 
such as the discontinuity on addition of an extra par- 
ticle to the system, which is reflected in the band gap 
problem |3| and the failure to correctly describe the re- 
sponse of a macroscopic system to an external electric 
field Q. This non-analytic behaviour is missing from 
the standard approximations to the Kohn-Sham poten- 
tial, i.e. the local density (LDA) and generalized gradient 
approximations (GGA), which may explain some of their 
limitations when applied to complex systems such as cat- 
alyzed chemical reactions ||. 

One way forward is to avoid the need to describe the 
non-analyticities by incorporating the true non-local na- 
ture of the exchange and correlation. Recently, a new 
realization of DFT, the generalized Kohn-Sham scheme 
(GKS) |||7]], has been proposed, in which the electrons 
move in an effective non-local potential. Among the GKS 
approaches, the screened-exchange LDA (sX-LDA) || 
appears to give the best performance, describing struc- 
tural properties with the same accuracy as the LDA but 
improving on its description of quasiparticle energies. 
However, this scheme is numerically as expensive as a 
standard Hartree-Fock calculation. 

In this work we propose a new generalized Kohn-Sham 
scheme, the E-GKS scheme, in which the exchange- 
correlation potential is intended to mimic the self- 
energy operator, which in many-body perturbation the- 
ory (MBPT) exactly describes exchange and correlation 
effects and, being truly non-local, is expected to be more 
amenable to approximation than the Kohn-Sham poten- 
tial. 

In the framework of MBPT, the total energy of a sys- 
tem of electrons moving in an external potential V ex t can 
be calculated by means of the Galitskii-Migdal formula: 
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where i and j are the indices of a suitable complete set of 
one-particle wavefunctions, h%j are the matrix elements 
of h = — V 2 /2 + V ex t and A is the spectral function of 
the Green's function. Atomic units are used throughout 
this Letter. 

^From a particular approximation for the self-energy, 
a model spectral function A can be obtained and the to- 
tal energy can be calculated via Eq. |. Using the GW 
approximation j|] , in its self-consistent formulation, total 
energies derived from Eq. [l] have recently been shown to 
be in very good agreement with quantum Monte Carlo 
data, when applied to 3D ||[lO| and 2D homogeneous 
electron gases. However this approximation is computa- 
tionally too expensive to be in competition with methods 
such as the LDA and GGA in applications to very large 
systems. The self-energy model that we propose will al- 
low us to calculate total energies within the same order 
of numerical cost as the LDA. 

We start by modelling the self-energy of a homoge- 
neous electron gas (jellium) of density no- For occupied 
states the energy- dependence of the self-energy is rela- 
tively weak B , suggesting the approximation [fy] : 



E(|r-r'|;e i )«£(|r-r / |;|i) 
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where /x is the chemical potential of the system. The 
self-energy for jellium at fi can be well reproduced by the 
following model function 



E°(r,r';no) = f(n )g(\r-r'\;n ) 
where f(n ) = V x L c DA (n ) and 

e — a(no)|r— r' 



g(|r-r'|;n ) = C(n )- 



(3) 



(4) 



The constant C(no) ensures that the chemical potential 
given by S° is correct and l/a(no) represents the range 
of the self-energy for a jellium of density no- 

In Fig. ^ the model function given by Eq. ^ is com- 
pared to the GW data of Ref. [Q , showing an excellent 
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agreement in the region in which E has the largest con- 
tribution. The small discrepancy that appears in the 
long-range region is unimportant for the total energy, 
owing to the oscillatory behaviour of the wavefunctions 
p3[ . In addition to its suitability for jellium, the function 
g defined in Eq. ^ has been shown to have the correct 
asymptotic decay in semiconducting systems Jl4[ . 

The key idea of this approach is to use a jellium-like 
self-energy to describe the exchange and correlation en- 
ergy of inhomogeneous systems. This is suggested by the 
fact that for several semiconductors the self-energy has 
been shown jl5) to be almost spherical and to have the 
same range as the self-energy of a jellium with the corre- 
sponding average density. A natural extension of Eq. ^ 
to an inhomogeneous system is: 

^(r,rMn]) = /( " (r)) + /(n(r/)) g (|r- r -|;n ) (5) 

where here no is the average density. This may be likened 
to a LDA exchange-correlation potential that has been 
made non-local by including a spreading function g. 

The self-energy in Eq. || is the self-energy of a system 
of non-interacting electrons, since it is energy indepen- 
dent and real. The method that we propose is equivalent 
to a mapping between the interacting system and the fic- 
titious non- interacting system represented by £°. If we 
choose the basis set <f>i to be formed by the wavefunc- 
tions of the quasiparticles in this fictitious sytem, then 
the spectral function A° is simply given by 



AUu>) = 6 ij 5(u-e°) 



(6) 



where e° are the quasiparticle energies of the non- 
interacting system. 

Using this spectral function A° and the Galitskii- 
Migdal formula (Eq. |]) the total energy of the model 
system can be obtained: 
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where, as usual, T is the kinetic energy, Eh is the Hartree 
energy, E ext is the external energy and 
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J dr #(r) J dr' E°(r, r', [n])^(r'). 



(8) 



As in other generalized Kohn-Sham schemes, a correc- 
tion term E ss is added so that the total energy given 
by the model is exact in the homogeneous limit. This 
term can be interpreted in MBPT as the remaining en- 
ergy associated with the satellite structure of the spectral 
function [fL6|. 
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where ei is the energy per particle obtained from Eq. 
in the limit of homogeneous densities: 



ei(n) = 
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where fcp(n) is the Fermi vector of a jellium of density n. 
The total energy of the interacting system is thus: 



E = T + E 
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This functional is minimized with respect to variations 
in the one-particle wavefunctions, yielding the following 
effective hamiltonian: 



1. 



h = --V' + V ext + V H 
where V ss (n(r)) — 5E SS /Sn(r) and 
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V m (r) = \f{n{r)) J dr' Rep(r, v')g{\v - r'|), 



(12) 
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in which p is the density matrix of the fictitious non- 
interacting system: 
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p(r,r') = > >*(r)<Mr') 



(14) 



In the limiting case of the spreading function g tending 
to a delta- function, the term i£° + V m + V ss reduces to 
the conventional LDA exchange-correlation potential, so 
that the LDA can be considered as a particular case of 
the more general scheme proposed here. 

The computational effort involved in solving the varia- 
tional equations determined by Eq. 12 is not significantly 
larger than in traditional DFT's. In a plane wave basis 
set, the matrix elements of the non-local potential take 
the simple form: 



(k + G|E°|k + G'> 

= /(G _ G , ) g (|k + G|)+ g (|k + G'| 



(15) 



i.e. they are just the product of the LDA exchange- 
correlation potential and an analytical function of the 
moduli of the vectors k+ G and k+ G'. The calculation 
of E° scales, as with local potentials, linearly with the 
number of reciprocal-lattice vectors ng. 

The potential V m defined by Eq. [l3] depends on 
the density matrix of the system rather than the den- 
sity itself, but it can be efficiently calculated in recip- 
rocal space, where the action of the non-local operator 
g on the wavefunctions is simply given by the prod- 
ucts g(G)(f>i(G). The calculation of V m also scales lin- 
early with no- In contrast, in Hartree- Fock or sX- 
LDA schemes, calculating the exchange-correlation ma- 
trix scales with the cube of tig- 
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By construction, the S-GKS is exact for homogeneous 
densities. In order to test its performance for inhomoge- 
neous systems, we have used it to calculate the linear re- 
sponse of a homogeneous electron gas to an external per- 
turbation 25V ex tCos(q ■ r). This is a stringent test since 
it involves the calculation of total energies for a whole 
family of systems with different values of the perturba- 
tion amplitude 5V ex t and the wavevector q. The total 
energy changes with respect to its unperturbed value by 
an amount x(q)^V e ^, t , where xil) is the response func- 
tion of the system, related to the local field factor G(q) 
via: 

/ s xo(q) / 1fi N 

X[q> l- W (q)[l-G(q)] X o(q) 1 ' 

where w(q) = An/q 2 is the Coulomb potential and Xo(q) 
is the non-interating Lindhard response function. For a 
given value of q, the total energy is calculated for sev- 
eral values of 5 V ex t and then the set of points {SV 2 xt , E) 
are fitted to a polynomial, whose first-order coefficient 
is x(q)- I n Fig. ||, the local field factor G(q) calculated 
within the S-GKS scheme is plotted against the wavevec- 
tor of the perturbation and compared to the quantum 
Monte Carlo (QMC) results of Moroni et al Q, the 
LDA and the GGA. For small values of q all methods 
agree well with QMC. The large q response, however, 
is not correctly described with the LDA, nor with the 
best GGA available, due to Perde, Burke and Ernzerhof 
(PBE) [Q. The E-GKS shows a significant improvement 
in the local field factor for large values of q. 

As a further test, the E-GKS scheme was applied to 
calculate structural properties of a typical semiconductor 
such as silicon. The bulk total energy was calculated us- 
ing an LDA pseudopotential jjlJ] due to Kerker and 
a plane wave basis set, and compared with the QMC re- 
sults of Ref . |EM , calculated with exactly the same pseu- 
dopotential. The results are shown in Table |. Both the 
S-GKS and the convential DFT's (LDA and GGA) total 
energies are in good agreement with the QMC results, al- 
though E-GKS agrees slightly less well [^5[ . The lattice 
parameter and bulk modulus are equally well described 
in the LDA, the GGA and the S-GKS scheme. 

Although describing quasiparticle spectra is not the 
primary goal of our scheme, we also calculated the quasi- 
particle energies of Si. The overall tendency is to cor- 
rectly increase the eigenvalues at the conduction band 
with respect to the LDA values. The direct band gap 
at r is increased from 2.6 eV in the LDA to 3.0 eV in 
our scheme, whereas the GW and experimental values 
are 3.4 eV. Valence-band widths in the E-GKS are in 
general overestimated. 

In summary, we have proposed a new method for cal- 
culating total energies at the same time as quasiparticle 
energies which is just as efficient as a standard DFT cal- 
culation, but that is constructed as a model of the self- 
energy, thus describing the exchange and correlation in 



a more realistic way than the LDA and GGA. This type 
of approach holds the prospect of enhanced accuracy of 
total-energy calculations by avoiding the pathological as- 
pects of the traditional Kohn-Sham exchange-correlation 
energy functional. 
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TABLE I. Total energy of bulk Si within the E-GKS scheme, the LDA and GGA (PW91) [22], calculated with a 4x4x4 
Monkhorst-Pack k-grid, and QMC [21], all with a cutoff of 15 Ry for the G- vectors. A LDA pseudopotential due to Kerker 
[20] was used in all four cases, which allows a comparison of the exchange-correlation functional under this particular choice 
of external potential. The lattice parameter and bulk modulus were also calculated, using a cutoff of 22 Ry for the G-points 
and a Monkhorst-Pack set of 14 special k-points in the Brillouin zone. The use of an LDA pseudopotential does not justify 
the comparison of our results with experiment, but nevertheless we present also the experimental values, for completeness. In 
common with the LDA and PW91, the E-GKS scheme shows generally good agreement with QMC and experiment. Note that 
quantities shown in parentheses have not been calculated with the pseudopotential of Ref. [20] . 



E (eV/atom) 



(rA) 



B (MBar) 



E-GKS 
LDA 
PW91 
QMC 
Expt. 



-107.65 
-107.90 
-108.17 
-108.3 b 



5.39 
5.39 
(5.46 a ) 
(5.45 c ) 
(5.43 d ) 



1.008 
0.967 
(0.92 a ) 
(1.03 c ) 
(0.992 d ) 



a From Ref. ffl 

b From Ref. 2U; estimated finite-size convergence error ±0.05 eV/atom 
c From Ref. g] 
d Cited in Ref. El 
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FIG. 1. Comparison between the model self-energy and the GW data from Ref. [12] for jellium of r, 



= 2. 
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FIG. 2. Local field factor for the linear response of jellium at r s = 2 in the E-GKS scheme compared to the QMC results 
[16], the LDA and the GGA of Ref. [17] (PBE) (which coincide). The numerical uncertainties in the E-GKS results are of the 
order of 0.05 and are not shown for clarity. The new scheme proposed considerably improves the local field factor for large 
values of q with respect to the LDA and to PBE. 
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